sep_smooths_year2
Load all required libraries.
library(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.3 v dplyr 1.0.7
## v tidyr 1.1.3 v stringr 1.4.0
## v readr 2.0.0 v forcats 0.5.1## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(plotly)##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutlibrary(broom)Read in raw data from RDS.
raw_data <- readRDS("./year2.RDS")Make a few small modifications to names and data for visualizations.
final_data <- raw_data %>% mutate(log_copy_per_L = log10(mean_copy_num_L)) %>%
rename(Facility = wrf) %>%
mutate(Facility = recode(Facility,
"NO" = "WRF A",
"MI" = "WRF B",
"CC" = "WRF C"))Seperate the data by gene target to ease layering in the final plot
#make three data layers
only_positives <<- subset(final_data, (!is.na(final_data$Facility)))
only_n1 <- subset(only_positives, target == "N1")
only_n2 <- subset(only_positives, target == "N2")
only_background <<-final_data %>%
select(c(date, cases_cum_clarke, new_cases_clarke, X7_day_ave_clarke)) %>%
group_by(date) %>% summarise_if(is.numeric, mean)
#specify fun colors
background_color <- "#7570B3"
seven_day_ave_color <- "#E6AB02"
marker_colors <- c("N1" = '#1B9E77',"N2" ='#D95F02')#remove facilty C for now
#only_n1 <- only_n1[!(only_n1$Facility == "WRF C"),]
#only_n2 <- only_n2[!(only_n2$Facility == "WRF C"),]
only_n1 <- only_n1[!(only_n1$Facility == "WRF A" & only_n1$date == "2020-11-02"), ]
only_n2 <- only_n2[!(only_n2$Facility == "WRF A" & only_n2$date == "2020-11-02"), ]Build the main plot
#first layer is the background epidemic curve
p1 <- only_background %>%
plotly::plot_ly() %>%
plotly::add_trace(x = ~date, y = ~new_cases_clarke,
type = "bar",
hoverinfo = "text",
text = ~paste('</br> Date: ', date,
'</br> Daily Cases: ', new_cases_clarke),
alpha = 0.5,
name = "Daily Reported Cases",
color = background_color,
colors = background_color,
showlegend = FALSE) %>%
layout(yaxis = list(title = "Clarke County Daily Cases", showline=TRUE)) %>%
layout(legend = list(orientation = "h", x = 0.2, y = -0.3))
#renders the main plot layer two as seven day moving average
p1 <- p1 %>% plotly::add_trace(x = ~date, y = ~X7_day_ave_clarke,
type = "scatter",
mode = "lines",
hoverinfo = "text",
text = ~paste('</br> Date: ', date,
'</br> Seven-Day Moving Average: ', X7_day_ave_clarke),
name = "Seven Day Moving Average Athens",
line = list(color = seven_day_ave_color),
showlegend = FALSE)
#renders the main plot layer three as positive target hits
p2 <- plotly::plot_ly() %>%
plotly::add_trace(x = ~date, y = ~mean_copy_num_L,
type = "scatter",
mode = "markers",
hoverinfo = "text",
text = ~paste('</br> Date: ', date,
'</br> Facility: ', Facility,
'</br> Target: ', target,
'</br> Copies/L: ', round(mean_copy_num_L, digits = 2)),
data = only_n1,
symbol = ~Facility,
marker = list(color = '#1B9E77', size = 8, opacity = 0.65),
showlegend = FALSE) %>%
plotly::add_trace(x = ~date, y = ~mean_copy_num_L,
type = "scatter",
mode = "markers",
hoverinfo = "text",
text = ~paste('</br> Date: ', date,
'</br> Facility: ', Facility,
'</br> Target: ', target,
'</br> Copies/L: ', round(mean_copy_num_L, digits = 2)),
data = only_n2,
symbol = ~Facility,
marker = list(color = '#D95F02', size = 8, opacity = 0.65),
showlegend = FALSE) %>%
layout(yaxis = list(title = "SARS CoV-2 Copies/L",
showline = TRUE,
type = "log",
dtick = 1,
automargin = TRUE)) %>%
layout(legend = list(orientation = "h", x = 0.2, y = -0.3))
#adds the limit of detection dashed line
p2 <- p2 %>% plotly::add_segments(x = as.Date("2021-06-30"),
xend = ~max(date + 10),
y = 3571.429, yend = 3571.429,
opacity = 0.35,
line = list(color = "black", dash = "dash")) %>%
layout(annotations = list(x = as.Date("2021-06-30"), y = 3.8, xref = "x", yref = "y",
text = "Limit of Detection", showarrow = FALSE))
p1 p2Combine the two main plot pieces as a subplot
#seperate n1 and n2 frames by site
#n1
wrf_a_only_n1 <- subset(only_n1, Facility == "WRF A")
wrf_b_only_n1 <- subset(only_n1, Facility == "WRF B")
wrf_c_only_n1 <- subset(only_n1, Facility == "WRF C")
#n2
wrf_a_only_n2 <- subset(only_n2, Facility == "WRF A")
wrf_b_only_n2 <- subset(only_n2, Facility == "WRF B")
wrf_c_only_n2 <- subset(only_n2, Facility == "WRF C")
#rejoin the old data frames then seperate in to averages for each plant.
wrfa_both <- full_join(wrf_a_only_n1, wrf_a_only_n2)%>%
select(c(date, mean_total_copies)) %>%
group_by(date) %>%
summarize_if(is.numeric, mean) %>%
ungroup() %>%
mutate(log_total_copies_both = log10(mean_total_copies))## Joining, by = c("date", "new_cases_clarke", "cases_cum_clarke", "X7_day_ave_clarke", "Facility", "collection_num", "target", "mean_copy_num_uL_rxn", "mean_copy_num_L", "sd_L", "mean_total_copies", "sd_total_copies", "log_copy_per_L")wrfb_both <- full_join(wrf_b_only_n1, wrf_b_only_n2)%>%
select(c(date, mean_total_copies)) %>%
group_by(date) %>%
summarize_if(is.numeric, mean) %>%
ungroup() %>%
mutate(log_total_copies_both = log10(mean_total_copies))## Joining, by = c("date", "new_cases_clarke", "cases_cum_clarke", "X7_day_ave_clarke", "Facility", "collection_num", "target", "mean_copy_num_uL_rxn", "mean_copy_num_L", "sd_L", "mean_total_copies", "sd_total_copies", "log_copy_per_L")wrfc_both <- full_join(wrf_c_only_n1, wrf_c_only_n2)%>%
select(c(date, mean_total_copies)) %>%
group_by(date) %>%
summarize_if(is.numeric, mean) %>%
ungroup() %>%
mutate(log_total_copies_both = log10(mean_total_copies))## Joining, by = c("date", "new_cases_clarke", "cases_cum_clarke", "X7_day_ave_clarke", "Facility", "collection_num", "target", "mean_copy_num_uL_rxn", "mean_copy_num_L", "sd_L", "mean_total_copies", "sd_total_copies", "log_copy_per_L")#get max date
maxdate <- max(wrfa_both$date)
mindate <- min(wrfa_both$date)Build loess smoothing figures figures
This makes the individual plots
#**************************************WRF A PLOT**********************************************
#add trendlines
#extract data from geom_smooth
#both extract
# *********************************span 0.6***********************************
#*****************Must always update the n = TOTAL NUMBER OF DAYS*************************
extract_botha <- ggplot(wrfa_both, aes(x = date, y = log_total_copies_both)) +
stat_smooth(aes(outfit=fit_botha<<-..y..), method = "loess", color = '#1B9E77',
span = 0.3, n = 169)## Warning: Ignoring unknown aesthetics: outfit#look at the fits to align dates and total observations
#both
extract_botha## `geom_smooth()` using formula 'y ~ x'
fit_botha## [1] 11.54419 11.57957 11.61484 11.65005 11.68522 11.72038 11.75556 11.79068
## [9] 11.82567 11.86056 11.89537 11.93015 11.96490 12.00006 12.03582 12.07191
## [17] 12.10803 12.14390 12.17925 12.21377 12.25095 12.28790 12.32119 12.35387
## [25] 12.38585 12.41701 12.44727 12.47654 12.50472 12.53258 12.56066 12.58848
## [33] 12.61558 12.64148 12.66435 12.68367 12.70098 12.71780 12.73565 12.75607
## [41] 12.78056 12.81066 12.84071 12.86980 12.90214 12.93612 12.97019 13.00275
## [49] 13.03222 13.05702 13.08297 13.11395 13.14546 13.17302 13.19213 13.20481
## [57] 13.21577 13.22459 13.23081 13.23399 13.23370 13.22950 13.21077 13.18928
## [65] 13.17468 13.15529 13.13223 13.10666 13.07970 13.05250 13.02619 13.00087
## [73] 12.97534 12.94884 12.92059 12.88981 12.85350 12.81082 12.76424 12.71622
## [81] 12.66923 12.62574 12.58822 12.55030 12.51307 12.47984 12.44544 12.41069
## [89] 12.37641 12.34345 12.31261 12.28473 12.25583 12.22403 12.19332 12.16768
## [97] 12.15111 12.14689 12.15316 12.16621 12.18231 12.19775 12.20880 12.21176
## [105] 12.21245 12.21312 12.21103 12.20946 12.20821 12.20708 12.20590 12.20447
## [113] 12.20260 12.20082 12.19949 12.19822 12.19657 12.19414 12.18754 12.17535
## [121] 12.16015 12.14448 12.13092 12.12202 12.12036 12.11931 12.11931 12.12793
## [129] 12.14266 12.16147 12.18231 12.20316 12.22199 12.23677 12.25170 12.27049
## [137] 12.29035 12.30852 12.32220 12.33470 12.35009 12.36677 12.38312 12.39755
## [145] 12.40844 12.41419 12.41641 12.41758 12.41737 12.41543 12.41143 12.40504
## [153] 12.39592 12.38462 12.37176 12.35710 12.34044 12.32153 12.30016 12.27610
## [161] 12.24958 12.22089 12.18992 12.15670 12.12125 12.08359 12.04376 12.00177
## [169] 11.95766#assign fits to a vector
both_trenda <- fit_botha
#extract y min and max for each
limits_botha <- ggplot_build(extract_botha)$data## `geom_smooth()` using formula 'y ~ x'limits_botha <- as.data.frame(limits_botha)
both_ymina <- limits_botha$ymin
both_ymaxa <- limits_botha$ymax
#reassign dataframes (just to be safe)
work_botha <- wrfa_both
#fill in missing dates to smooth fits
work_botha <- work_botha %>% complete(date = seq(min(date), max(date), by = "1 day"))
date_vec_botha <- work_botha$date
#create a new smooth dataframe to layer
smooth_frame_botha <- data.frame(date_vec_botha, both_trenda, both_ymina, both_ymaxa)#WRF A
#plot smooth frames
p_wrf_a <- plotly::plot_ly() %>%
plotly::add_lines(x = ~date_vec_botha, y = ~both_trenda,
data = smooth_frame_botha,
hoverinfo = "text",
text = ~paste('</br> Date: ', date_vec_botha,
'</br> Median Log Copies: ', round(both_trenda, digits = 2)),
line = list(color = '#1B9E77', size = 8, opacity = 0.65),
showlegend = FALSE) %>%
layout(xaxis = list(range = c(mindate - 7, maxdate + 7))) %>% #buffer here
plotly::add_ribbons(x ~date_vec_botha, ymin = ~both_ymina, ymax = ~both_ymaxa,
showlegend = FALSE,
opacity = 0.25,
hoverinfo = "text",
text = ~paste('</br> Date: ', date_vec_botha, #leaving in case we want to change
'</br> Max Log Copies: ', round(both_ymaxa, digits = 2),
'</br> Min Log Copies: ', round(both_ymina, digits = 2)),
name = "",
fillcolor = '#1B9E77',
line = list(color = '#1B9E77')) %>%
layout(yaxis = list(title = "Total Log10 SARS CoV-2 Copies",
showline = TRUE,
automargin = TRUE)) %>%
layout(xaxis = list(title = "Date")) %>%
layout(title = "WRF A") %>%
plotly::add_markers(x = ~date, y = ~log_total_copies_both,
data = wrfa_both,
hoverinfo = "text",
showlegend = FALSE,
text = ~paste('</br> Date: ', date,
'</br> Actual Log Copies: ', round(log_total_copies_both, digits = 2)),
marker = list(color = '#1B9E77', size = 6, opacity = 0.65))
p_wrf_asave(p_wrf_a, file = "./site_objects/wrf_a_year2.rda")#**************************************WRF B PLOT**********************************************
#add trendlines
#extract data from geom_smooth
#both extract
# *********************************span 0.6***********************************
#*****************Must always update the n = TOTAL NUMBER OF DAYS*************************
extract_bothb <- ggplot(wrfb_both, aes(x = date, y = log_total_copies_both)) +
stat_smooth(aes(outfit=fit_bothb<<-..y..), method = "loess", color = '#D95F02',
span = 0.3, n = 169)## Warning: Ignoring unknown aesthetics: outfit#look at the fits to align dates and total observations
#both
extract_bothb## `geom_smooth()` using formula 'y ~ x'
fit_bothb## [1] 10.70850 10.80013 10.88992 10.97793 11.06425 11.14897 11.23215 11.31393
## [9] 11.39426 11.47298 11.54993 11.62495 11.69787 11.76903 11.83872 11.90677
## [17] 11.97299 12.03721 12.09923 12.15888 12.21597 12.27113 12.32672 12.38354
## [25] 12.44028 12.49563 12.54828 12.59692 12.64024 12.68208 12.72528 12.76705
## [33] 12.80461 12.83517 12.85445 12.86248 12.86294 12.85951 12.85589 12.85576
## [41] 12.86280 12.86763 12.87170 12.88259 12.89429 12.90635 12.91833 12.92979
## [49] 12.94028 12.94936 12.96105 12.97732 12.99466 13.00953 13.01842 13.02261
## [57] 13.02560 13.02720 13.02721 13.02544 13.02170 13.01580 12.99953 12.98221
## [65] 12.96853 12.94781 12.92302 12.89714 12.87315 12.85403 12.84275 12.83973
## [73] 12.84160 12.84560 12.84900 12.84903 12.85034 12.85759 12.86808 12.87912
## [81] 12.88802 12.89209 12.88863 12.89043 12.89043 12.87727 12.86137 12.84292
## [89] 12.82211 12.79912 12.77415 12.74738 12.71011 12.65926 12.60352 12.55154
## [97] 12.51201 12.47984 12.44525 12.40987 12.37531 12.34322 12.31520 12.29290
## [105] 12.29283 12.29328 12.28281 12.28360 12.29036 12.29779 12.30058 12.29344
## [113] 12.27106 12.23095 12.17883 12.12169 12.06657 12.02047 11.96964 11.90127
## [121] 11.82382 11.74574 11.67547 11.62147 11.59219 11.54492 11.50357 11.50717
## [129] 11.52124 11.54379 11.57290 11.60659 11.64292 11.67992 11.73901 11.82972
## [137] 11.93225 12.02683 12.09366 12.15116 12.22579 12.30989 12.39576 12.47573
## [145] 12.54213 12.58726 12.61048 12.61942 12.61892 12.61378 12.60881 12.60884
## [153] 12.61868 12.63257 12.64229 12.64897 12.65372 12.65765 12.66188 12.66752
## [161] 12.67211 12.67570 12.68039 12.68495 12.68924 12.69312 12.69647 12.69913
## [169] 12.70099#assign fits to a vector
both_trendb <- fit_bothb
#extract y min and max for each
limits_bothb <- ggplot_build(extract_bothb)$data## `geom_smooth()` using formula 'y ~ x'limits_bothb <- as.data.frame(limits_bothb)
both_yminb <- limits_bothb$ymin
both_ymaxb <- limits_bothb$ymax
#reassign dataframes (just to be safe)
work_bothb <- wrfb_both
#fill in missing dates to smooth fits
work_bothb <- work_bothb %>% complete(date = seq(min(date), max(date), by = "1 day"))
date_vec_bothb <- work_bothb$date
#create a new smooth dataframe to layer
smooth_frame_bothb <- data.frame(date_vec_bothb, both_trendb, both_yminb, both_ymaxb)#WRF B
#plot smooth frames
p_wrf_b <- plotly::plot_ly() %>%
plotly::add_lines(x = ~date_vec_bothb, y = ~both_trendb,
data = smooth_frame_bothb,
hoverinfo = "text",
text = ~paste('</br> Date: ', date_vec_bothb,
'</br> Median Log Copies: ', round(both_trendb, digits = 2)),
line = list(color = '#D95F02', size = 8, opacity = 0.65),
showlegend = FALSE) %>%
layout(xaxis = list(range = c(mindate - 7, maxdate + 7))) %>% #buffer here
plotly::add_ribbons(x ~date_vec_bothb, ymin = ~both_yminb, ymax = ~both_ymaxb,
showlegend = FALSE,
opacity = 0.25,
hoverinfo = "text",
text = ~paste('</br> Date: ', date_vec_bothb, #leaving in case we want to change
'</br> Max Log Copies: ', round(both_ymaxb, digits = 2),
'</br> Min Log Copies: ', round(both_yminb, digits = 2)),
name = "",
fillcolor = '#D95F02',
line = list(color = '#D95F02')) %>%
layout(yaxis = list(title = "Total Log10 SARS CoV-2 Copies",
showline = TRUE,
automargin = TRUE)) %>%
layout(xaxis = list(title = "Date")) %>%
layout(title = "WRF B") %>%
plotly::add_markers(x = ~date, y = ~log_total_copies_both,
data = wrfb_both,
hoverinfo = "text",
showlegend = FALSE,
text = ~paste('</br> Date: ', date,
'</br> Actual Log Copies: ', round(log_total_copies_both, digits = 2)),
marker = list(color = '#D95F02', size = 6, opacity = 0.65))
p_wrf_bsave(p_wrf_b, file = "./site_objects/wrf_b_year2.rda")#**************************************WRF C PLOT********************************************** #add trendlines #extract data from geom_smooth # *********************************span 0.6*********************************** #*****************Must always update the n = TOTAL NUMBER OF DAYS*************************
extract_bothc <- ggplot(wrfc_both, aes(x = date, y = log_total_copies_both)) +
stat_smooth(aes(outfit=fit_bothc<<-..y..), method = "loess", color = '#E7298A',
span = 0.3, n = 169)## Warning: Ignoring unknown aesthetics: outfit#look at the fits to align dates and total observations
#both
extract_bothc## `geom_smooth()` using formula 'y ~ x'
fit_bothc## [1] 10.69069 10.78286 10.87218 10.95842 11.04133 11.12068 11.19623 11.26809
## [9] 11.33665 11.40206 11.46447 11.52404 11.58092 11.63315 11.67949 11.72121
## [17] 11.75958 11.79585 11.83129 11.86717 11.88980 11.91011 11.94208 11.97408
## [25] 12.00586 12.03718 12.06781 12.09750 12.12600 12.15082 12.17105 12.18881
## [33] 12.20624 12.22549 12.24782 12.27254 12.29882 12.32582 12.35269 12.37860
## [41] 12.40271 12.42887 12.45463 12.47425 12.48965 12.50202 12.51252 12.52233
## [49] 12.53261 12.54454 12.55605 12.56518 12.57325 12.58158 12.59149 12.60791
## [57] 12.63246 12.66155 12.69156 12.71888 12.73990 12.75101 12.75737 12.76093
## [65] 12.75539 12.74162 12.72198 12.69885 12.67458 12.65155 12.63212 12.61233
## [73] 12.58829 12.56254 12.53756 12.51588 12.49386 12.46744 12.43851 12.40893
## [81] 12.38060 12.35540 12.33519 12.32252 12.31049 12.30055 12.29965 12.30384
## [89] 12.30919 12.31173 12.30752 12.29262 12.26362 12.22425 12.18127 12.14139
## [97] 12.11136 12.08332 12.04746 12.00758 11.96748 11.93096 11.90183 11.88388
## [105] 11.87170 11.86184 11.86063 11.86808 11.88147 11.89807 11.91516 11.93001
## [113] 11.93989 11.95748 11.99045 12.02853 12.06150 12.07909 12.08318 12.08305
## [121] 12.07959 12.07368 12.06620 12.05805 12.05009 12.03022 12.00849 11.98994
## [129] 11.96004 11.92359 11.88538 11.85021 11.82286 11.80812 11.80308 11.80105
## [137] 11.80123 11.80277 11.80487 11.81483 11.83705 11.86634 11.89751 11.92536
## [145] 11.94471 11.95036 11.94212 11.92546 11.90373 11.88030 11.85852 11.84173
## [153] 11.83330 11.83023 11.82734 11.82468 11.82232 11.82030 11.81869 11.81753
## [161] 11.81676 11.81622 11.81594 11.81599 11.81640 11.81721 11.81843 11.82009
## [169] 11.82222#assign fits to a vector
both_trendc <- fit_bothc
#extract y min and max for each
limits_bothc <- ggplot_build(extract_bothc)$data## `geom_smooth()` using formula 'y ~ x'limits_bothc <- as.data.frame(limits_bothc)
both_yminc <- limits_bothc$ymin
both_ymaxc <- limits_bothc$ymax
#reassign dataframes (just to be safe)
work_bothc <- wrfc_both
#fill in missing dates to smooth fits
work_bothc <- work_bothc %>% complete(date = seq(min(date), max(date), by = "1 day"))
date_vec_bothc <- work_bothc$date
#create a new smooth dataframe to layer
smooth_frame_bothc <- data.frame(date_vec_bothc, both_trendc, both_yminc, both_ymaxc)#WRF C
#plot smooth frames
p_wrf_c <- plotly::plot_ly() %>%
plotly::add_lines(x = ~date_vec_bothc, y = ~both_trendc,
data = smooth_frame_bothc,
hoverinfo = "text",
text = ~paste('</br> Date: ', date_vec_bothc,
'</br> Median Log Copies: ', round(both_trendc, digits = 2)),
line = list(color = '#E7298A', size = 8, opacity = 0.65),
showlegend = FALSE) %>%
layout(xaxis = list(range = c(mindate - 7, maxdate + 7))) %>% #buffer here
plotly::add_ribbons(x ~date_vec_bothc, ymin = ~both_yminc, ymax = ~both_ymaxc,
showlegend = FALSE,
opacity = 0.25,
hoverinfo = "text",
text = ~paste('</br> Date: ', date_vec_bothc, #leaving in case we want to change
'</br> Max Log Copies: ', round(both_ymaxc, digits = 2),
'</br> Min Log Copies: ', round(both_yminc, digits = 2)),
name = "",
fillcolor = '#E7298A',
line = list(color = '#E7298A')) %>%
layout(yaxis = list(title = "Total Log10 SARS CoV-2 Copies",
showline = TRUE,
automargin = TRUE)) %>%
layout(xaxis = list(title = "Date")) %>%
layout(title = "WRF C") %>%
plotly::add_markers(x = ~date, y = ~log_total_copies_both,
data = wrfc_both,
hoverinfo = "text",
showlegend = FALSE,
text = ~paste('</br> Date: ', date,
'</br> Actual Log Copies: ', round(log_total_copies_both, digits = 2)),
marker = list(color = '#E7298A', size = 6, opacity = 0.65))
p_wrf_csave(p_wrf_c, file = "./site_objects/wrf_c_year2.rda")keeping in case
#save(wrfa_both, file = "./plotly_objs/wrfa_both.rda")
#save(wrfb_both, file = "./plotly_objs/wrfb_both.rda")
#save(wrfc_both, file = "./plotly_objs/wrfc_both.rda")#save(date_vec_botha, file = "./plotly_objs/date_vec_botha.rda")
#save(date_vec_bothb, file = "./plotly_objs/date_vec_bothb.rda")
#save(date_vec_bothc, file = "./plotly_objs/date_vec_bothc.rda")#save(both_ymina, file = "./plotly_objs/both_ymina.rda")
#save(both_ymaxa, file = "./plotly_objs/both_ymaxa.rda")
#save(both_yminb, file = "./plotly_objs/both_yminb.rda")
#save(both_ymaxb, file = "./plotly_objs/both_ymaxb.rda")
#save(both_yminc, file = "./plotly_objs/both_yminc.rda")
#save(both_ymaxc, file = "./plotly_objs/both_ymaxc.rda")